We follow chapter II of the book and model the simplest case, adding features as we go along.
In [1]:
# We define an abstract type AgentInfo with two subtypes: NoAgent and Agent.
# This allows for removing an agent from the board while keeping his entry in the agent list.
# Also it allows for recursive type definition, an unresolved issue in Julia #269
abstract AgentInfo
type NoAgent <: AgentInfo
end
In [2]:
type Place
x::Int
y::Int
agent::AgentInfo
sugar::Int
capacity::Int
end
In [3]:
type Agent <: AgentInfo
place::Place
vision::Int
metabolism::Int
stash::Int
end
In [4]:
# To keep record of the simulation, we keep all agent steps in memory for now
# Alternatively, only the current agent list could be kept or all agents with
# the lattice at each point.
type AgentList
agents::Vector{AgentInfo}
end
type Scape
lattice::Array{Place}
steps::Vector{AgentList}
end
In [5]:
# convenience functions
agents(Scape) = last(Scape.steps).agents
Out[5]:
In [6]:
# We define sugarscape capacity similar to graphs in the book
function init_capacity()
grid = 50
hump1 = [15, 15]
hump2 = [40, 40]
dist(xy, hump) = hypot(xy[1]-hump[1], xy[2]-hump[2])
capacity = Int[max(0, 4 - ifloor(min(dist([i,j], hump1), dist([i,j], hump2))/5)) for i=1:grid,j=1:grid]
end
Out[6]:
In [7]:
using PyPlot
In [8]:
#plot capacity (not yet in lattice)
capacity = init_capacity()
cm = ColorMap([Color.RGB(1,1,1), Color.RGB(1,1,0)],5,1.5)
ind_i, ind_j = collect(ind2sub(size(capacity),1:length(capacity)))
scatter(ind_i, ind_j,s=capacity*20, c=capacity,cmap=cm,marker="o",lw=0)
xlim(0, 50); ylim(0, 50)
title("capacity")
Out[8]:
In [9]:
using StatsBase: sample
In [10]:
function init_scape(capacity; N_agents=400)
gridx, gridy = size(capacity)
# create empty lattice with sugar levels at full capacity
lattice = [Place(i,j, NoAgent(), capacity[i,j], capacity[i,j]) for i=1:gridx, j=1:gridy]
#populate with agents
agents = AgentInfo[]
for place in sample(lattice, N_agents, replace=false)
vision = rand(1:6)
metabolism = rand(1:4)
stash = 10 #not specified in book
agent = Agent(place, vision, metabolism, stash)
push!(agents, agent)
place.agent = agent
end
Scape(lattice, [AgentList(agents)])
end
Out[10]:
In [12]:
@time scape = init_scape(init_capacity());
In [13]:
# convenience functions to retrieve only alive agents
alive(agent::AgentInfo) = isa(agent, Agent)
function living(scape)
allagents = agents(scape)
allagents[map(alive, allagents)]
end
Out[13]:
Let's also define scape plotting
In [14]:
function plot(scape::Scape)
lattice = scape.lattice
gridx, gridy = size(lattice)
sugar = reshape([p.sugar for p in lattice], gridx, gridy)
s_i, s_j = collect(ind2sub(size(lattice), 1:length(lattice)))
colmap = ColorMap([Color.RGB(1,1,1), Color.RGB(1,1,0)],5,1.5)
PyPlot.scatter(s_i, s_j, s=20*sugar, c=sugar, cmap=cm, marker="o", lw=0)
a_i = [a.place.x for a in living(scape)]
a_j = [a.place.y for a in living(scape)]
scatter(a_i, a_j, c="red", lw=0)
gridx
xlim(0, gridx+1); ylim(0, gridy+1)
title("Sugarscape step $(length(scape.steps))")
end
Out[14]:
In [15]:
plot(scape)
Out[15]:
In [16]:
pygui(true)
Out[16]:
In [17]:
@time plot(scape)
Out[17]:
An agent can jump as far as she can see according to her vision, but only in the four main directions.
In [18]:
#functions to select neighbouring place with circular boundaries
left(p::Place, lat::Array{Place}) = lat[p.x==1 ? end : p.x-1, p.y]
right(p::Place, lat::Array{Place}) = lat[p.x==end ? 1 : p.x+1, p.y]
up(p::Place, lat::Array{Place}) = lat[p.x, p.y==1 ? end : p.y-1]
down(p::Place, lat::Array{Place}) = lat[p.x, p.y==end ? 1 : p.y+1]
Out[18]:
In [19]:
# convenience function
Base.isempty(place::Place) = isa(place.agent, NoAgent)
Out[19]:
In [20]:
function move(agent::Agent, dest_place::Place)
agent.place === dest_place && return
isempty(dest_place) || error("destination place occupied")
agent.place.agent = NoAgent()
dest_place.agent = agent
agent.place = dest_place
return nothing
end
Out[20]:
In [42]:
function view_place(agent::Agent, step_function::Function, lattice::Array{Place})
p = agent.place
p_out = pstep = p #if no unoccupied place in sight, stay put
sugar = -1 #Agent will move to empty place rather than stay put
for _ in 1:agent.vision
pstep = step_function(pstep, lattice)
# Only move further if more sugar, not on equal sugar level.
# As in book, agents will take steps of size=1 when far from sugar
# mountains, but this could be improved to larger steps.
# It seems from Animation II-1 that agents do not jump (as a horse), but
# only slide like a rook.
if isempty(pstep)
if pstep.sugar > sugar
sugar = pstep.sugar
p_out = pstep
end
else #disallow jumps
break
end
end
p_out
end
Out[42]:
In [22]:
# to sort array of places by sugar, define `isless`
Base.isless(p1::Place, p2::Place) = p1.sugar < p2.sugar
Out[22]:
In [23]:
function harvest(agent)
agent.place.sugar == 0 && return(nothing)
agent.stash += agent.place.sugar
agent.place.sugar = 0
nothing
end
Out[23]:
In [43]:
function move(scape::Scape)
lat = scape.lattice
@inbounds for agent in shuffle(living(scape))
view_places = Place[]
# Randomize search directions in case of equal sugar place, as in book.
for wind in shuffle([left, right, up, down])
push!(view_places, view_place(agent, wind, lat))
end
select_place = last(sort!(view_places))
move(agent, select_place)
harvest(agent)
end
end
Out[43]:
In [44]:
# test methods
move(scape)
PyPlot.clf()
plot(scape)
Out[44]:
In [70]:
function grow(scape::Scape; α=1)
@inbounds for place in scape.lattice
place.sugar = min(place.capacity, place.sugar + α)
end
end
Out[70]:
In [71]:
grow(scape)
PyPlot.clf()
plot(scape)
Out[71]:
In [72]:
function consume(scape::Scape)
@inbounds for agent in living(scape)
agent.stash -= agent.metabolism
end
end
Out[72]:
In [73]:
function Base.kill(scape::Scape)
a = agents(scape)
for i = 1:length(a)
if alive(a[i]) && a[i].stash < 0
a[i].place.agent = NoAgent()
a[i] = NoAgent()
end
end
end
Out[73]:
In [74]:
consume(scape)
In [75]:
function timestep(scape::Scape)
grow(scape)
move(scape)#includes harvestig
consume(scape)
kill(scape)
push!(scape.steps, AgentList(agents(scape)))
nothing
end
Out[75]:
In [90]:
pygui(true)
PyPlot.clf()
scape = init_scape(init_capacity())
plot(scape)
Out[90]:
In [95]:
@time timestep(scape)
PyPlot.clf()
plot(scape)
Out[95]:
In [96]:
@time for cnt = 1:100
timestep(scape)
if mod(cnt, 10) == 0
sleep(.01)
PyPlot.clf() #clear figure from artifacts
plot(scape)
end
end
In [ ]:
In [82]:
function run(n)
scape = init_scape(init_capacity())
for cnt = 1:n
timestep(scape)
end
scape
end
Out[82]:
In [83]:
run(1);
In [85]:
@time scape = run(1000);
In [86]:
length(scape.steps)
Out[86]: